Regresia segmentată și cauzalitatea

De la observație la inferență cauzală

Claudiu Papasteri

Prezentarea

Urmăriți prezentarea pe:

https://quarto.org/docs/presentations/

sau scanând codul


Un pattern fără context (1)

click pentru a vedea codul
# Packages and settings
packs <- c("ggplot2", "papaja")
success <- suppressWarnings(sapply(packs, require, character.only = TRUE))
install.packages(names(success)[!success])
sapply(names(success)[!success], require, character.only = TRUE)

theme_set(papaja::theme_apa())  # set APA style theme for plots

# Get data from package
if(!require(MultiKink)){install.packages("MultiKink")}; library(MultiKink)
data <- MultiKink::triceps
data <- subset(data, age <= 30)

# Plot the data
plot <- ggplot(data, aes(x = age, y = triceps)) +
  geom_point(alpha = .7, fill = "black", color = "grey50") + 
  xlab("Vârstă") +
  ylab("TSF")

Investigați codul:

  • Pachete & Setări
  • Datele provin dintr-un pachet
  • Se generează un grafic

Un pattern fără context (2)

Ce vedeți?

Cum am putea modela?

click pentru a vedea codul
plot

Un pattern fără context (3)

Ce reprezintă linia? Este pattern-ul descris bine de o relație liniară?

click pentru a vedea codul
# Linear Regression
# plot + geom_smooth(formula = y ~ x, method = "lm", se = TRUE, color = "red") # fast way to plot it
lm_lin <- lm(triceps ~ age, data = data)   # fit regression model and then plot it
pred_lin <- predict(lm_lin, se.fit = TRUE, data = data) 

plot_lin <- 
  plot +
  geom_line(aes(y = pred_lin$fit), size = 1, color = "red") +
  geom_ribbon(aes(ymax = pred_lin$fit + pred_lin$se.fit, ymin = pred_lin$fit - pred_lin$se.fit), 
              fill = "red", alpha = .3) +
  geom_text(x = 5, y = 35, color = "red", 
            label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_lin)$r.squared, 3))))

suppressWarnings(print(plot_lin))

Un pattern fără context (4)

Este pattern-ul mai bine descris de o relație curbilinie (polinomială)? Care dintre trenduri, pătratic (roșu) sau cubic (albastru) se potrivește mai bine datelor?

click pentru a vedea codul
# Polynomial Regression
lm_quad <- lm(triceps ~ age + I(age^2), data = data)
lm_cub <- lm(triceps ~ age + I(age^2) + I(age^3), data = data)

pred_quad <- predict(lm_quad, se.fit = TRUE, data = data) 
pred_cub <- predict(lm_cub, se.fit = TRUE, data = data) 

plot_poly <-
  plot +
  geom_line(aes(y = pred_quad$fit), size = 1, color = "blue") +
  geom_ribbon(aes(ymax = pred_quad$fit + pred_quad$se.fit, ymin = pred_quad$fit - pred_quad$se.fit), 
              fill = "blue", alpha = .1) +
  geom_text(x = 5, y = 35, color = "blue", 
            label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_quad)$r.squared, 3)))) +
  geom_line(aes(y = pred_cub$fit), size = 1, color = "red") +
  geom_ribbon(aes(ymax = pred_cub$fit + pred_cub$se.fit, ymin = pred_cub$fit - pred_cub$se.fit), 
              fill = "red", alpha = .1) +
  geom_text(x = 5, y = 30, color = "red", 
            label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_cub)$r.squared, 3))))

suppressWarnings(print(plot_poly))

Un pattern fără context (5)

Reinventăm regresia segmentată!

  1. Împărțim datele în grupe de vârstă (segmentăm variabila \(x\)).

  2. Potrivim regresii liniare pe fiecare set de date rezultat.

click pentru a vedea codul
plot

Un pattern fără context (6)

Reinventăm regresia segmentată!

\[ y = \beta_{0} + \beta_{1}x + \epsilon \;\;\;\; pentru \;\;\;\; x \leq k_{1} \\ y = \beta_{0} + \beta_{1}x + \epsilon \;\;\;\; pentru \;\;\;\; k_{1} < x \leq k_{2} \\ ... \\ y = \beta_{0} + \beta_{1}x + \epsilon \;\;\;\; pentru \;\;\;\; x \geq k_{p} \\ \]

Denumim \(k\) punctul de întrerupere (nod).

Rezultă \(p\) segmente.

Un pattern fără context (7)

Reinventăm regresia segmentată, în cod!

click pentru a vedea codul
# Piecewise regression
data_k1 <- data[data$age < 5, ]
data_k2 <- data[data$age >= 5 & data$age < 10, ]
data_k3 <- data[data$age >= 10 & data$age < 20, ]   
data_k4 <- data[data$age >= 20, ]
  
lm_p1 <- lm(triceps ~ age, data = data_k1)
lm_p2 <- lm(triceps ~ age, data = data_k2)
lm_p3 <- lm(triceps ~ age, data = data_k3)
lm_p4 <- lm(triceps ~ age, data = data_k4)

pred_p1 <- predict(lm_p1, se.fit = TRUE, data = data_k1)
pred_p2 <- predict(lm_p2, se.fit = TRUE, data = data_k2)
pred_p3 <- predict(lm_p3, se.fit = TRUE, data = data_k3)
pred_p4 <- predict(lm_p4, se.fit = TRUE, data = data_k4)

plot_pwr <- 
  plot + 
  geom_line(data = data_k1, aes(y = pred_p1$fit), size = 1, color = "red") +
  geom_ribbon(data = data_k1, aes(ymax = pred_p1$fit + pred_p1$se.fit, ymin = pred_p1$fit - pred_p1$se.fit),
              fill = "red", alpha = .1) + 
  geom_line(data = data_k2, aes(y = pred_p2$fit), size = 1, color = "red") +
  geom_ribbon(data = data_k2, aes(ymax = pred_p2$fit + pred_p2$se.fit, ymin = pred_p2$fit - pred_p2$se.fit),
              fill = "red", alpha = .1) +
  geom_line(data = data_k3, aes(y = pred_p3$fit), size = 1, color = "red") +
  geom_ribbon(data = data_k3, aes(ymax = pred_p3$fit + pred_p3$se.fit, ymin = pred_p3$fit - pred_p3$se.fit),
              fill = "red", alpha = .1) +
  geom_line(data = data_k4, aes(y = pred_p4$fit), size = 1, color = "red") +
  geom_ribbon(data = data_k4, aes(ymax = pred_p4$fit + pred_p4$se.fit, ymin = pred_p4$fit - pred_p4$se.fit),
              fill = "red", alpha = .1)  

suppressWarnings(print(plot_pwr))

Un pattern fără context (8)

Reinventăm regresia segmentată, în cod!

Ce observați despre segmentele din Regresia Segmentată?

Un pattern fără context (9)

De fapt, ceea ce este denumit în literatură Regresie Segmentată nu este Regresie Segmentată, ci Regresie cu Spline-uri Liniare.

Spline-urile permit o interpolare netedă și continuă între punctele de întrerupere (noduri). Între noduri este calculată o regresie (polinomială).

Un pattern fără context (10)

Re-reinventăm regresia segmentată!

  1. Împărțim datele în grupe de vârstă (segmentăm variabila \(x\)) => \(k\) noduri.

  2. Folosim o funcție treaptă încât să obținem o singură ecuație de regresie care păstrează doar un intercept \(\beta_{0}\) și câte o pantă pentru fiecare set de date rezultat.

\[ y = \beta_{0} + \beta_{1}x + \beta_{2}(x - k_{1}) + \beta_{3}(x - k_{2}) + ... + \beta_{p-1}(x - k_{p}) + \epsilon \\ \begin{equation} Unde \;\;\;\; (x - k)_{+} = \begin{cases} 0 \;\;\;\;\;\;\;\;\; dacă \; x < k \\ x - k \;\;\; dacă \; x \geq k \\ \end{cases} \end{equation} \]

Denumim \(k\) punctul de întrerupere (nod).

Rezultă \(p\) segmente.

Un pattern fără context (11)

Re-reinventăm regresia segmentată, în cod!

click pentru a vedea codul
# Linear splines
# Add the terms (x − k) when x ≥ k . 
# We will do this by adding I ((age − k)∗(age >= k)) terms to the linear model. 
# Note that (age >= k) is a logical statement that will be 0 (FALSE) or 1 (TRUE) and I() evaluates that all expression.
lm_linsp <- lm(triceps ~ age + 
                 I((age - 5) * (age >= 5)) +
                 I((age - 10) * (age >= 10)) +
                 I((age - 20) * (age >= 20)),
               data = data)
pred_linsp <- predict(lm_linsp, se.fit = TRUE, data = data)

plot_linsp <-                  
  plot +
  geom_line(data = data, aes(y = pred_linsp$fit, x = age), size = 1, col = "red") +
  geom_ribbon(data = data, aes(ymax = pred_linsp$fit + pred_linsp$se.fit, ymin = pred_linsp$fit - pred_linsp$se.fit), 
              fill = "red", alpha = .1) +
  geom_text(x = 5, y = 35, color = "red", 
            label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_linsp)$r.squared, 3))))

suppressWarnings(print(plot_linsp))

Un pattern fără context (12)

Bonus: Mai sus scria că regresiile între noduri sunt polinomiale, iar polinomul poate, desigur, să nu fie de ordinul 1.

Optimizând numărul \(k\) de noduri și ordinul polinomului obținem o metodă foarte robustă, deși nesofisticată, de învățare a pattern-urilor din date. Un exemplu de Regresie cu Spline-uri Polinomiale de ordin 2:

click pentru a vedea codul
# Quadratic splines
lm_quadsp <- lm(triceps ~ age + 
                  I((age - 5) * (age >= 5)) + I((age - 5)^2 * (age >= 5)) + 
                  I((age - 10) * (age >= 10)) + I((age - 10)^2 * (age >= 10)) + 
                  I((age - 20) * (age >= 20)) + I((age - 20)^2 * (age >= 20)), 
                data = data)

pred_quadsp <- predict(lm_quadsp, se.fit = TRUE, data = data)

plot_quadsp <-
  plot +
  geom_line(data = data, aes(y = pred_quadsp$fit, x = age), size = 1, col = "red") +
  geom_ribbon(data = data, aes(ymax = pred_quadsp$fit + pred_quadsp$se.fit, ymin = pred_quadsp$fit - pred_quadsp$se.fit), 
              fill = "red", alpha = .1) +
  geom_text(x = 5, y = 35, color = "red", 
            label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_quadsp)$r.squared, 3))))

suppressWarnings(print(plot_quadsp))

Contextul datelor, elucidarea misterului

Măsurăm constructe:

  • TSF = grosimea pliului cutanat al tricepsului, o metrică economică și convenabilă pentru evaluarea obezității.

    • Validitate: față de alți indici are avantajul de a reprezenta distribuţia grăsimii.

    • Utilitate predictivă: la fiecare creștere cu 1 mm a TSF, riscul de deces scade cu 4%, riscul de deces din cauze cardiovasculare cu 6% (Li et al., 2022).

Vârsta codează timpul în unități discrete dacă am presupune că indivizii sunt interșanjabili. Totuși, multe variabile se asociază cu TSF (ex. genul, nutriția etc), iar imaginea transversală nu spune povestea filmului longitudinal.

Date longitudinale în Psihologie

Deoarece cea mai mare parte a psihologiei implică așa-numitele procese non-ergodice, cercetarea pe eșantione mari nu poate oferi informații fidele despre procese la nivel individual. (Hamaker, 2012).

  • O descriere simplificată a non-ergotismului: individul, în timp, nu obține rezultatul mediu al grupului.

Monitorizare în timp real a sistemelor idiografice

  • Pentru o discuție a utilizării conceptului și complianței cu metodologia în context psihoterapeutic vedeți (Schiepek et al., 2016).

Datele longitudinale

Datele disponibile de la Journal of Open Psychology Data: (Heino, 2022)

  • Monitorizare dinamică cognitivă (test Stroop) după mindfulness
  • Un singur individ peste 900 de zile, meditație zilnică de 20 minute

Repozitoriu OSF
Datele se găsesc și în Github-ul prezentării.

Lupta cu date longitudinale (1)

click pentru a vedea codul
# Packages and settings
packs <- c("ggplot2", "papaja", "osfr", "tidyverse", "scales", "lspline", "plotly")
success <- suppressWarnings(sapply(packs, require, character.only = TRUE))
install.packages(names(success)[!success])
sapply(names(success)[!success], require, character.only = TRUE)

theme_set(papaja::theme_apa())  # set APA style theme for plots

# Read data from folder -- here we will read directly from OSF repository
# csv_name <- "cognitive_dynamics_heino.csv"
# csv_path <- here::here("datasets", csv_name) 
# df <- read.csv(csv_path)   

# Read data from OSF repository
repo_url <- "https://osf.io/w9v28/"
csv_name <- "cognitive_dynamics_heino.csv"
repo <- osfr::osf_retrieve_node(repo_url); Sys.sleep(10)
repo_files <- osfr::osf_ls_files(repo)
# we download to a temporary file that gets deleted after
temp_dir <- tempdir()
osfr::osf_download(repo_files[repo_files$name == csv_name, ], path = temp_dir); Sys.sleep(10)
csv_path <- file.path(temp_dir, csv_name)
df <- read.csv(csv_path)
unlink(csv_path)   

# Transform
df_clean <- 
  df %>%
  # keep only days with record of date
  tidyr::drop_na(date) %>%
  # convert to ISO 8601 date format, format used is not consistent
  tidyr::separate(date, into = c("day", "month", "year"), sep = "\\.", remove = FALSE) %>%
  tidyr::separate(time_of_day, into = c("hour", "minute"), sep = "\\:", remove = FALSE) %>%
  # there are some typos, exclude them
  dplyr::mutate(across(c(day, month, year, hour, minute), as.numeric)) %>%
  dplyr::filter(year >= 2014, year <= 2017, month < 12, day < 31) %>%
  dplyr::mutate(
    date = lubridate::ymd(paste(year, month, day, sep = ' ')),
    date_posixct = as.POSIXct(date),
    hm = lubridate::hm(paste(hour, minute, sep = ' '))
  ) %>%
  dplyr::filter(date > lubridate::as_date("2014-07-01")) %>%
  dplyr::select(rowNumber, date, date_posixct, hm, stroop_incongruent_post_meditation) %>%
  tidyr::drop_na(stroop_incongruent_post_meditation) %>%
  dplyr::mutate(y = stroop_incongruent_post_meditation,
                moments = dplyr::row_number())

# Plot
plot <- 
  ggplot(df_clean, aes(x = date_posixct, y = y)) +
  geom_point() +
  scale_x_datetime(date_labels = "%Y-%m", breaks= scales::date_breaks("6 month")) +
  xlab("") +
  ylab("Timp răspuns Stroop inc. Post mindfulness") + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

Investigați codul:

  • Pachete & Setări
  • Datele sunt citite din repozitorul OSF
  • Datele sunt curățate și transformate
  • Se generează un grafic

Lupta cu date longitudinale (2)

Ce observați?

Regresia Liniară

Putem modela mai bine?

click pentru a vedea codul
# Linear Regression
lm_lin <- lm(y ~ date_posixct, data = df_clean)
pred_lin <- predict(lm_lin, se.fit = TRUE, data = df_clean)

plot_lin <-                  
  plot +
  geom_line(data = df_clean, aes(y = pred_lin$fit, x = date_posixct), size = 1, col = "blue") +
  geom_ribbon(data = df_clean, aes(ymax = pred_lin$fit + pred_lin$se.fit, ymin = pred_lin$fit - pred_lin$se.fit), 
              fill = "blue", alpha = .1) +
  geom_text(x = as.POSIXct("2014-11-30"), y = 21, color = "blue", 
            label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_lin)$r.squared, 3))))

suppressWarnings(print(plot_lin))

Regresia segmentată

De această dată nu o vom mai coda “manual” (folosind doar base R), ci vom folosi pachetul lspline pentru a obține noduri în urma împărțirii datelor în intervale egale de timp.

click pentru a vedea codul
# Linear Splines
# Instead of doing it manually, we will use `lspline` package
# knots set manually lspline()
# knots set at breaks dividing the range of x into q equal-frequency intervals qlspline()
# knots set at breaks dividing the range of x into n equal-width intervals elspline()
lm_linsp2 <- lm(y ~ lspline::qlspline(date_posixct, q = 2), data = df_clean)
lm_linsp3 <- lm(y ~ lspline::qlspline(date_posixct, q = 3), data = df_clean)
lm_linsp4 <- lm(y ~ lspline::qlspline(date_posixct, q = 4), data = df_clean)
lm_linsp5 <- lm(y ~ lspline::qlspline(date_posixct, q = 5), data = df_clean)
lm_linsp6 <- lm(y ~ lspline::qlspline(date_posixct, q = 6), data = df_clean)
lm_linsp10 <- lm(y ~ lspline::qlspline(date_posixct, q = 10), data = df_clean)

df_clean$knots2 <- predict(lm_linsp2, se.fit = TRUE, data = df_clean)$fit
df_clean$knots3 <- predict(lm_linsp3, se.fit = TRUE, data = df_clean)$fit
df_clean$knots4 <- predict(lm_linsp4, se.fit = TRUE, data = df_clean)$fit
df_clean$knots5 <- predict(lm_linsp5, se.fit = TRUE, data = df_clean)$fit
df_clean$knots6 <- predict(lm_linsp6, se.fit = TRUE, data = df_clean)$fit
df_clean$knots10 <- predict(lm_linsp10, se.fit = TRUE, data = df_clean)$fit

Interacționați cu datele

Câte noduri să folosim?

click pentru a vedea codul
plot_interact <- 
  plot_ly(data = df_clean, x = ~date_posixct) %>% 
  add_markers(y = ~y, showlegend = FALSE, name = "") %>% 
  add_lines(x = ~date_posixct, y = ~knots2, visible = "legendonly", name = "2 noduri") %>%
  add_lines(x = ~date_posixct, y = ~knots3, visible = "legendonly", name = "3 noduri") %>%
  add_lines(x = ~date_posixct, y = ~knots4, visible = "legendonly", name = "4 noduri") %>%
  add_lines(x = ~date_posixct, y = ~knots5, visible = "legendonly", name = "5 noduri") %>%
  add_lines(x = ~date_posixct, y = ~knots6, visible = "legendonly", name = "6 noduri") %>%
  layout(
    xaxis = list(
      title = "", 
      tickformat = "%Y-%m"
  ))

plot_interact

Câte noduri?

click pentru a vedea codul
# Plot PWR with different number of knots
plot_knots <- 
  plot +
  geom_line(data = df_clean, aes(y = knots5, x = date_posixct), alpha = .8, size = 1, color = "red1") +
  geom_text(x = as.POSIXct("2014-11-30"), y = 21, color = "red1", 
          label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_linsp5)$r.squared, 3)))) + 
  geom_line(data = df_clean, aes(y = knots6, x = date_posixct), alpha = .8, size = 1, color = "red3") +
  geom_text(x = as.POSIXct("2014-11-30"), y = 19.5, color = "red3", 
            label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_linsp6)$r.squared, 3)))) +
  geom_line(data = df_clean, aes(y = knots10, x = date_posixct), alpha = .8, size = 1, color = "blue") +
  geom_text(x = as.POSIXct("2014-11-30"), y = 17.5, color = "blue", 
            label = latex2exp::TeX(paste0("$\\R^2 = $", round(summary(lm_linsp10)$r.squared, 3))))

suppressWarnings(print(plot_knots))

Pe ce bazăm inferența

Inferență bazată pe model

VS.

Inferență bazată pe design

Cauzalitatea (după Hume)

Cauza și efectul sunt interconectate

Cauza precede (temporal) efectul

Cauza și efectul covariază consistent

Nu există alte explicații alternative

Regresia segmentată și cauzalitatea (1)

Design-ul de Discontinuitate în Regresie (RDD)

  • cvasi-experimental pre-post
  • manipularea variabilei independente poate apărea natural în mediu => inferență cauzală din date observaționale
  • NU vom discuta despre ea la acest curs, dar puteți citi mai multe

Regresia segmentată și cauzalitatea (2)

Studiile Experimentale cu un Singur Caz (SCED)

  • manipularea variabilei independente
  • măsurători repetate
  • calitatea dovezilor cauzale echivalentă cu cea din studii clinice randomizate (RCT)

Regresia segmentată în SCED

Istoria modelului

      

Istoria modelului pentru SCED

 
Model ANOVA  

 
Model cu o pantă globală  

 
Model regresie segmentată  

Istoria modelului în cod (1)

Codați cu mine pașii:

  • Pachete & Setări
  • Date
  • Calcul statistici pentru grafice
click pentru a vedea codul
# Packages and settings
packs <- c("tidyverse", "papaja", "latex2exp", "cowplot")
success <- suppressWarnings(sapply(packs, require, character.only = TRUE))
install.packages(names(success)[!success])
sapply(names(success)[!success], require, character.only = TRUE)

if(!require(ggbrace)){devtools::install_github("nicolash2/ggbrace")}; library(ggbrace)

theme_set(papaja::theme_apa())  # set APA style theme for plots


# Data
df <- tibble(
  y = c(c(3, 4, 6, 5, 4), c(8, 7, 8, 9, 9)),
  AB = c(rep("A", 5), rep("B", 5)), 
  t = 1:10
)

# Stats
df_A <- df[df$AB == "A", ]
df_B <- df[df$AB == "B", ]  
mean_A <- mean(df_A$y, na.rm = TRUE)
mean_B <- mean(df_B$y, na.rm = TRUE)
intercept_A <- summary(lm(y ~ t, data = df_A))$coefficients[1, 1]
slope_A <- summary(lm(y ~ t, data = df_A))$coefficients[2, 1]
intercept_B <- summary(lm(y ~ t, data = df_B))$coefficients[1, 1]
slope_B <- summary(lm(y ~ t, data = df_B))$coefficients[2, 1]

Istoria modelului în cod (2)

Plotarea datelor, un schelet pentru restul graficelor.

click pentru a vedea codul
# Plot Skeleton
plot_skeleton <- 
  ggplot() +
  geom_point(data = df, aes(x = t, y = y)) +
  geom_rect(aes(xmin = 5.4, xmax = 5.6, ymin = 1, ymax = 10), 
            fill = "black", alpha = 0.2) +
  geom_line(data = df_A, aes(x = t, y = y)) +
  geom_line(data = df_B, aes(x = t, y = y)) +
  scale_x_continuous(breaks = 1:10, limits = c(1, 10)) +
  scale_y_continuous(breaks = 1:10, limits = c(1, 10)) +
  coord_equal()

plot_skeleton

Istoria modelului în cod (3)

Mediile celor două faze.

click pentru a vedea codul
# Plot Means
plot_skeleton +    
  geom_segment(aes(x = 1, xend = 5.5, y = mean_A, yend = mean_A), 
               color = "tan1") +
  geom_segment(aes(x = 5.5, xend = 10, y = mean_B, yend = mean_B), 
               color = "tan1")

Istoria modelului în cod (4)

Panta globală (peste datele din ambele faze).

click pentru a vedea codul
# Plot Global Slope
plot_skeleton + 
  stat_smooth(data = df, aes(x = t, y = y), geom = "line",
              formula = y ~ x, method = "lm", se = FALSE, 
              linetype = "longdash", size = 1, alpha = 0.6, 
              color = "yellow4")

Istoria modelului în cod (5)

Pantele individuale ale celor două faze.

click pentru a vedea codul
# Plot Phase Slopes
plot_skeleton +
  stat_smooth(data = df_A, aes(x = t, y = y), geom = "line", 
              formula = y ~ x, method = "lm", se = FALSE,
              alpha = 0.7, size = 1, color = "yellow3") +
  stat_smooth(data = df_B, aes(x = t, y = y), geom = "line",
              formula = y ~ x, method = "lm", se = FALSE,
              alpha = 0.7, size = 1, color = "yellow3")

Istoria modelului în cod (6)

Ce arată graficul?

Istoria modelului în cod (7)

Cum am obținut graficul?

# Plot all together
plot_skeleton +
  geom_segment(aes(x = 1, xend = 5.5, y = mean_A, yend = mean_A), 
               color = "tan1") +
  geom_segment(aes(x = 5.5, xend = 10, y = mean_B, yend = mean_B), 
               color = "tan1") +  
  stat_smooth(data = df, aes(x = t, y = y), geom = "line",
              formula = y ~ x, method = "lm", se = FALSE, 
              linetype = "longdash", size = 1, alpha = 0.6, 
              color = "yellow4") +
  stat_smooth(data = df_A, aes(x = t, y = y), geom = "line", 
              formula = y ~ x, method = "lm", se = FALSE,
              alpha = 0.7, size = 1, color = "yellow3") +
  stat_smooth(data = df_B, aes(x = t, y = y), geom = "line",
              formula = y ~ x, method = "lm", se = FALSE,
              alpha = 0.7, size = 1, color = "yellow3")

Modelul ANOVA

\(y = \beta_{0} + \beta_{1} D + \epsilon\)

Modelul ANOVA

  • diferență în medii (nivel) = \(\beta_{1}\)

Modelul cu o Pantă Globală

\(y = \beta_{0} + \beta_{1} D + \beta_{2} t + \epsilon\)

Modelul cu o Pantă Globală

  • (Gorsuch, 1983)
  • diferență în medii (nivel) = \(\beta_{1}\), dar doar dacă nu există trend (adică \(\beta_{2}\) = 0)
  • ia în calcul trendul, dar presupune trend identic între faze

Modelul de Regresie Segmentată

\(y = \beta_{0} + \beta_{1} D + \beta_{2} t + \beta_{3} D(t-k) + \epsilon\)

Modelul de Regresie Segmentată

  • (Center et al., 1985)
  • diferență în medii (nivel) = \(\beta_{1}D\)
  • trend baseline = \(\beta_{2}t\)
  • termen de interacțiune ce cuantifică diferența în pantă între cele două faze = \(\beta_{3} D(t-k)\)
  • unde k = nr. observații în faza A (baseline)

Interpretare

\[ y = \color{blue} { \beta_{0} } + \color{green} {\beta_{1} D } + \color{blue} { \beta_{2} t } + \color{lime} { \beta_{3} D(t-k) } + \epsilon \]

\(\color{blue} { \beta_{0} }\) este interceptul (scorul evaluat la \(t = 0\) și \(D = 0\)), deci indică valoarea de început a liniei de regresie în faza A.

\(\color{green} {\beta_{1} }\) poate fi interpretat ca schimbarea nivelului dintre faza A și faza B necounfounduit cu posibilele efecte de trend, mai precis, este diferența dintre scorurile prezise ale ambelor linii de regresie la prima măsurătoare a fazei B.

\(\color{blue} { \beta_{2} }\) captează trendul din faza A.

\(\color{lime} { \beta_{3} }\), care este de obicei parametrul cel mai de interes, reprezintă modificarea în trend de la faza A la faza B.

Datele SCED

Datele provin din celebrul studiu realizat de (Singh et al., 2007), digitalizate de Rumen Manolov.

Variabile:

tier Id numeric participant
id Nume participant
time Index al momentului măsurătorii (\(t_{0}\) = 1).
phase Variabilă dummy: 0 tratamentul nu a început (baseline), 1 tratamentul a început.
score_physical Scorul participantului la agresivitate fizică.
score_verbal Scorul participantului la agresivitate verbală.


Download
Repozitoriu OSF
Datele participantului “Michael” utilizate în această prezentare se găsesc și în Github-ul prezentării.

Preambul în cod

click pentru a vedea codul
# Packages and settings
packs <- c("ggplot2", "papaja", "scda", "scan")
success <- suppressWarnings(sapply(packs, require, character.only = TRUE))
install.packages(names(success)[!success])
sapply(names(success)[!success], require, character.only = TRUE)

theme_set(papaja::theme_apa())  # set APA style theme for plots

# Get data
# Michael's data from `scda` package (can be downloaded from presentation Gihub)
data <- scda::Singh
data_Michael <- subset(data, id == "Michael")

# Plot
plot <- 
  ggplot() +
    geom_point(data = data_Michael, aes(x = time, y = score_physical)) +
    geom_rect(aes(xmin = 4.4, xmax = 4.6, ymin = 0, ymax = 3), 
              fill = "black", alpha = 0.2) +
    geom_line(data = data_Michael[data_Michael$phase == 0, ], aes(x = time, y = score_physical)) +
    geom_line(data = data_Michael[data_Michael$phase == 1, ], aes(x = time, y = score_physical)) +
    scale_x_continuous(breaks = 1:18, limits = c(1, 18)) +
    scale_y_continuous(breaks = 0:3, limits = c(0, 3)) +
    xlab("Timp") +
    ylab("Scor fizic") +
    coord_equal()

Investigați codul:

  • Pachete & Setări
  • Datele provin dintr-un pachet
  • Se generează un grafic

Vizualizarea datelor

Ce observați?

click pentru a vedea codul
plot

Analize cu pachetul scda

click pentru a vedea codul
# Package scda -----------------------------------------------------------
# Piecewise regression
pwr_scda <- scda::piecewiseRegr(data = data_Michael,
                           timeVar = "time",
                           yVar = "score_physical",
                           phaseVar = "phase",
                           robust = FALSE,
                           outputWidth = 1.2,
                           outputHeight = 1.2,
                           yRange = c(0,7),
                           showPlot = FALSE)

Grafic cu pachetul scda

click pentru a vedea codul
pwr_scda$output$plot

Rezultate cu pachetul scda

click pentru a vedea codul
pwr_scda
Piecewise Regression Analysis (N = 18)

Model statistics:

  Model deviance:                         6.41
  R squared for null model:               0.63
  R squared for test model:               0.7
  R squared based effect size:            0.19
  Effect size (delta_t):                  3.36
  Standardized effect size (delta_ts):    5.47
  Effect size (delta):                    0.74
  Standardized effect size (delta_s):     1.2

  Effect size evaluated at point:         18

Regression coefficients:

  Intercept:        [1.49; 3.91] (point estimate = 2.7)
  Level change:    [-2.42; 1.42] (point estimate = -0.5)
  Trend phase 1:   [-0.95; 0.35] (point estimate = -0.3)
  Change in trend: [-0.47; 0.85] (point estimate = 0.19)

Mărime a efectului cu pachetul scda

\(delta_{t}\) a lui (Swaminathan et al., 2014) compară scorurile prezise la un anumit punct de măsurare (t) din faza B cu predicția din faza A extrapolată la t.

click pentru a vedea codul
# Effect size
# Swaminathan et al. (2014) compares the predicted scores from the PWR model at a given measurement point (t) in phase B with the prediction from the A phase extrapolated to t.
pwr_scda$output$delta_t 
(Intercept) 
   3.361538 

Regresie segmentată generalizată cu pachetul scda

Extindere a modelului utilizată în design-uri mai complexe (ex. ABA cu retragere, ABAB cu inversiune).

click pentru a vedea codul
# Generalized PWR 
pwrgen_scda <- scda::genPwr(data = data_Michael,
                           timeVar = "time",
                           yVar = "score_physical",
                           phaseVar = "phase")

pwrgen_scda
pwrgen_scda$output$plot

Analize cu pachetul scan

click pentru a vedea codul
# Package scan -----------------------------------------------------------
# Transform into scdf list format
# data_scdf <- as.list(split(data, f = as.factor(data$id)))
data_Michael_scdf <- scan::scdf(
  values = data_Michael$score_physical, 
  covariate = data_Michael$score_verbal,
  mt = data_Michael$time,
  phase = LETTERS[data_Michael$phase + 1],        # data_Michael$phase
  name = "Michael"
)

# Piecewise regression
# "B&L-B" starts the dummy var at 1, "H-M" starts at 0 which is also used by scda::piecewiseRegr
pwr_scan <- scan::plm(data = data_Michael_scdf, 
                      dvar = "values", pvar = "phase", mvar = "mt", 
                      model = "H-M",                                    
                      AR = 0, family = "gaussian",
                      trend = TRUE, level = TRUE, slope = TRUE, r_squared = TRUE)

Grafic cu pachetul scan

click pentru a vedea codul
plot(data_Michael_scdf)

Rezultate cu pachetul scan

click pentru a vedea codul
pwr_scan
Piecewise Regression Analysis

Dummy model:  H-M 

Fitted a gaussian distribution.
F(3, 14) = 11.07; p = 0.001; R² = 0.703; Adjusted R² = 0.640

                  B   2.5% 97.5%    SE      t     p delta R²
Intercept      3.00  1.376 4.624 0.829  3.620 0.003         
Trend mt      -0.30 -0.893 0.293 0.303 -0.991 0.338   0.0208
Level phase B -0.50 -2.258 1.258 0.897 -0.557 0.586   0.0066
Slope phase B  0.19 -0.409 0.790 0.306  0.621 0.544   0.0082

Autocorrelations of the residuals
 lag    cr
   1  0.32
   2 -0.12
   3 -0.24

Formula: values ~ 1 + mt + phaseB + interB

Autocorelație și covariate cu pachetul scan

click pentru a vedea codul
# Autocorrelation
pwr_scan_autocor <- scan::autocorr(data = data_Michael_scdf, 
  dvar = "values", pvar = "phase", mvar = "mt", lag_max = 3)

# Multivariate piecewise regression
mvpwr_scan <- scan::mplm(data = data_Michael_scdf, 
  dvar = c("values", "covariate"), pvar = "phase", mvar = "mt", 
  model = "H-M", trend = TRUE, level = TRUE, slope = TRUE)

mvpwr_scan
Multivariate piecewise linear model

Dummy model: H-M 

Coefficients: 
              values covariate
(Intercept)     3.00     6.500
Trend          -0.30    -1.000
Level Phase B  -0.50     0.271
Slope Phase B   0.19     0.870

Formula: y ~ 1 + mt + phaseB + interB

Type III MANOVA Tests: Pillai test statistic
              Df test stat approx F num Df den Df Pr(>F)   
(Intercept)    1     0.558     8.22      2     13 0.0049 **
Trend          1     0.175     1.38      2     13 0.2868   
Level Phase B  1     0.058     0.40      2     13 0.6782   
Slope Phase B  1     0.148     1.13      2     13 0.3524   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The following variables were used in this analysis:
'values/ covariate' as dependent variable, 'phase' as phase variable, and 'mt' as measurement-time variable.

Reproducem fără pachete

Totul în base R.

click pentru a vedea codul
# Do it in base R
k_time <- 5   # (t-k_time) starts with 0 in B phase => k is number of first time of B phase
pwr <- lm(score_physical ~ 1 + phase + time + I(phase * (time - k_time) * (time >= k_time)), data = data_Michael)
# here phase and (time >= k_time) are both 0 in A and 1 in B ... don't need both dummy vars
summary(pwr)$coefficients  # same results as scan

pred_pwr <- predict(pwr, se.fit = TRUE, data = data_Michael)

plot_base <- 
  ggplot(data = data_Michael, aes(x = time, y = score_physical)) +
    geom_point(alpha = 0.55, color = "black") + 
    geom_line(aes(y = pred_pwr$fit), size = 1, color = "red") +
    geom_ribbon(aes(ymax = pred_pwr$fit + pred_pwr$se.fit, ymin = pred_pwr$fit - pred_pwr$se.fit), fill = "red", alpha = .1)  

# scda centers t so that it begins with 0 in both phases => difference in intercept between scda and scan/base solutions
# to get same result for intercept as scda:
data_Michael2 <- data_Michael
data_Michael2$time <- data_Michael2$time - 1
k_time <- 4
pwr_likescda <- lm(score_physical ~ 1 + phase + time + I(phase * (time - k_time) * (time >= k_time)), data = data_Michael2)
summary(pwr_likescda)$coefficients  # same as scda

Reproducem fără pachete

Observăm o diferență în estimarea interceptului între scda și scan datorată centrării timpului. Puteți citi mai multe despre recomandări de codificare a timpului în SCED în (Huitema & Mckean, 2000) (mai ales, schema standard H-M ce poartă numele autorilor).

În base R putem reproduce ambele rezultate.

La fel ca scda

                                                Estimate Std. Error    t value
(Intercept)                                    2.7000000  0.5661223  4.7692871
phase                                         -0.5000000  0.8969429 -0.5574491
time                                          -0.3000000  0.3026051 -0.9913910
I(phase * (time - k_time) * (time >= k_time))  0.1901099  0.3059124  0.6214521
                                                  Pr(>|t|)
(Intercept)                                   0.0002993711
phase                                         0.5860238465
time                                          0.3383210448
I(phase * (time - k_time) * (time >= k_time)) 0.5442892970

La fel ca scan

                                                Estimate Std. Error    t value
(Intercept)                                    3.0000000  0.8287183  3.6200482
phase                                         -0.5000000  0.8969429 -0.5574491
time                                          -0.3000000  0.3026051 -0.9913910
I(phase * (time - k_time) * (time >= k_time))  0.1901099  0.3059124  0.6214521
                                                 Pr(>|t|)
(Intercept)                                   0.002785467
phase                                         0.586023846
time                                          0.338321045
I(phase * (time - k_time) * (time >= k_time)) 0.544289297

Analize în GUI

https://manolov.shinyapps.io/Regression/

Bibliografie

Center, B. A., Skiba, R. J., & Casey, A. (1985). A Methodology for the Quantitative Synthesis of Intra-Subject Design Research. The Journal of Special Education, 19(4), 387–400. https://doi.org/10.1177/002246698501900404
Gorsuch, R. L. (1983). Three methods for analyzing limited time-series (n of 1) data. Behavioral Assessment, 5, 141–154.
Hamaker, E. L. (2012). Why researchers should think "within-person": A paradigmatic rationale. In Handbook of research methods for studying daily life. (pp. 43–61). The Guilford Press.
Heino, M. T. J. (2022). Cognitive Dynamics of a Single Subject: 1428 Stroop Tests and Other Measures in a Mindfulness Meditation Context Over 2.5 Years. Journal of Open Psychology Data, 10. https://doi.org/10.5334/jopd.51
Huitema, B. E., & Mckean, J. W. (2000). Design Specification Issues in Time-Series Intervention Models. Educational and Psychological Measurement, 60(1), 38–58. https://doi.org/10.1177/00131640021970358
Li, W., Yin, H., Chen, Y., Liu, Q., Wang, Y., Qiu, D., Ma, H., & Geng, Q. (2022). Associations between adult triceps skinfold thickness and all-cause, cardiovascular and cerebrovascular mortality in NHANES 19992010: A retrospective national study. Frontiers in Cardiovascular Medicine, 9. https://doi.org/10.3389/fcvm.2022.858994
Schiepek, G., Aichhorn, W., Gruber, M., Strunk, G., Bachler, E., & Aas, B. (2016). Real-time monitoring of psychotherapeutic processes: Concept and compliance. Frontiers in Psychology, 7. https://doi.org/10.3389/fpsyg.2016.00604
Singh, N. N., Lancioni, G. E., Winton, A. S. W., Adkins, A. D., Wahler, R. G., Sabaawi, M., & Singh, J. (2007). Individuals with Mental Illness Can Control their Aggressive Behavior Through Mindfulness Training. Behavior Modification, 31(3), 313–328. https://doi.org/10.1177/0145445506293585
Swaminathan, H., Rogers, H. J., & Horner, R. H. (2014). An effect size measure and Bayesian analysis of single-case designs. Journal of School Psychology, 52(2), 213–230. https://doi.org/10.1016/j.jsp.2013.12.002